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We present new sequences of general relativistic, quasiequilibrium black hole-neutron star binaries. 
We solve for the gravitational field in the conformal thin-sandwich decomposition of Einstein's field 
equations, coupled to the equations of relativistic hydrostatic equilibrium for a perfect fluid. We 
account for the black hole by solving these equations in the background metric of a Schwarzschild 
black hole whose mass is much larger than that of the neutron star. The background metric is treated 
in Kerr-Schild as well as isotropic coordinates. For the neutron star, we assume a polytropic equation 
of state with adiabatic index F = 2, and solve for both irrotational and corotational configurations. 
By comparing the results of irrotational and synchronized configurations with the same background 
metric, we conclude that the effect of the rotation on the location of tidal break-up is only on the 
order of a few percent. The different choices in the background also lead to differences of order 
a few percent, which may be an indication of the level to which these configurations approximate 
quasiequilibrium. 

PACS numbers: 04.30.Db, 04.25. Dm, 04.40.Dg 



I. INTRODUCTION 

Black hole-neutron star (hereafter BHNS) binary sys- 
tems are, together with other compact binaries, among 
the most promising sources of gravitational waves for de- 
tection by ground-based interferometers such as LIGO 
[|, GEO600 @, TAMA300 0, and VIRGO [J, or for 
the planned space-based mission LISA [|. Since theo- 
retical predictions for the waveforms from such compact 
binaries are needed both for the identification and inter- 
pretation of any astrophysical signals, significant effort 
has gone into the theoretical modeling of these binaries 
and their inspiral (see [1,0] for recent reviews). 

The orbital separation of a BHNS binary decreases 
as energy and angular momentum are dissipated by the 
emission of gravitational radiation, which also has the ef- 
fect of circularizing the orbit, until the system eventually 
coalesces. This coalescence can take two qualitatively dif- 
ferent forms. If the neutron star reaches the innermost 
stable circular orbit (hereafter ISCO) while still stable 
against tidal disruption, it will likely plunge and fall into 
the black hole promptly. Alternately, if the neutron star 
is tidally disrupted outside the ISCO, the final fate may 
be a black hole surrounded by a disk. This model is 
a candidate central engine for short-period gamma-ray 
bursts, since the high efficiency of accretion onto a black 
hole can help to explain the huge luminosities seen in 
these sources [IUCkl]. Additionally, we may obtain in- 
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formation about the equation of state (hereafter EOS) 
of matter at nuclear densities through the detection of a 
gravitational wave signal, because the characteristic fre- 
quency at the tidal disruption separation is in the most 
sensitive range of the ground-based detectors [ll| . 

The difference in the final fate is determined primar- 
ily by the mass ratio of the black hole to the neutron 
star. When the black hole mass -Mbh is much larger 
than the neutron star mass A/ns, the tidal force from 
the black hole is not so large compared with the neutron 
star self-gravity to disrupt the neutron star, even at the 
ISCO. Indeed, the ISCO occurs at a radius proportional 
to the black hole mass, whereas the Roche limit separa- 
tion, which would mark the beginning of mass transfer, 
1/3 . 

scales like Af BH - in the limit of an extremely large mass 
ratio [Hjl- Thus, the larger the black hole mass, the fur- 
ther outward the ISCO lies relative to the Roche limit 
separation. On the other hand, a stellar mass black hole 
will have a tidal field sufficiently large with respect to 
the neutron star self-gravity to deform and disrupt the 
companion outside the ISCO. The critical value of the 
mass ratio falls in the range Mbh /A/ns ~ 4 for reason- 
able neutron star models but varies depending on 
assumptions about the neutron star EOS (see for a 
thorough discussion). 

Until now, much effort has been devoted to the com- 
putation of the BHNS binary systems. Most dynamical 
simulations so far have been carried out in a Newtonian 
framework [1 El, El, El El, El HI , but see [fl HJ for 
approximate relativistic treatments. Quasiequilibrium 
models of BHNS binaries have been constructed adopt- 
ing various different approximations. Several authors (in- 
cluding [12]) have modeled BHNS binaries as Newtonian 
ellipsoids around point masses, generalizing the classic 
Roche model for incompressible stars [23|. These ellip- 
soidal calculations have also been gene ralized to include 
relativistic effects [H, S3, HI [Mill- Recently, these 
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models have been generalized to include higher order de- 
formations than ellipsoid by expanding the background 
black hole metric to higher order than quadrupole [28j]. 
Equilibrium models of Newtonian BHNS binaries also 
have been constructed by solving the exact fluid equa- 
tions numerically and again treating the black hole as 
a point-mass (29[. With the possible exception of [30| . 
Baumgarte et al. [H[ (hereafter BSS) so far provide the 
only self-consistent relativistic treatment of BHNS bina- 
ries in quasiequilibrium. 

In BSS, BHNS binaries are constructed by solving the 
constraint equations of general relativity, decomposed in 
the conformal thin-sandwich formalism, together with 
the Euler equation for the neutron star matter, which 
takes an algebraic form if the system can be taken to 
be stationary in a corotating frame. The equations are 
solved in the background of a Schwarzschild black hole, 
which accounts for the neutron star's companion. 

In this paper we generalize the findings of BSS in two 
ways. As discussed above we allow of irrotational instead 
of corotational fluid flow, which is more realistic astro- 
physically [12, [Hj]. To do so we adopt the formalism 
for constructing irrotational stars as developed in [34l ]. 
We also generalize the results of BSS by expressing the 
black hole background in a different coordinate system. 
In BSS, this background was expressed in Kerr-Schild 
coordinates, which we compare here with a black hole 
background expressed in isotropic coordinates. Since the 
two coordinate system represent different slicings of the 
Schwarzschild space-time, the resulting initial data are 
physically distinct solutions of the constraint equations. 
Finally, in this paper we adopt a numerical code that is 
based on the spectral method numerical libraries known 
as LORENE [35j (as opposed to the finite difference im- 
plementation of BSS). We maintain several of the other 
assumptions made by BSS, including extreme mass ratio 
Mbh ^ -Mns and polytropic equations of state, and plan 
to relax these in the future. 

The spectral techniques used here have previously been 
employed to compute quasiequilibrium sequences of bi- 
nary neutron stars (H, [H, HE Hi HI S3 . There are sev- 
eral advantages in using spectral methods. One is that it 
is easy to treat spherical coordinates, which are suitable 
for solving figures like stars, and to fit the coordinates 
to trace the stellar surface. This fitting technique plays 
an important role in solving irrotational configurations. 
Another advantage is that it is possible to achieve more 
rapid convergence compared to finite difference methods, 
until configurations reach a separation very close to tidal 
break-up and the appearance of discontinuities in physi- 
cal quantities (For more details on the spectral methods 
techniques we use, we refer the readers to [H, l42j.) 
In our code, the set of equations we solve for the gravita- 
tional field is essentially equivalent to that of BSS, while 
the hydrostatic equations are the same as those found in 

The paper is organized as follows. In Section II, we 
briefly summarize the formulation and explain the so- 



lution procedure. The tests of the numerical code are 
shown in Section III, and the results are presented in Sec- 
tion IV. In Section V we discuss the effects of the choice 
of the black hole background metric, and we summarize 
in Section VI. 

Throughout the present paper, we adopt geometrical 
units, G = c = 1, where G denotes the gravitational 
constant and c the speed of light, respectively. Latin 
and Greek indices denote purely spatial and space-time 
components, respectively. 

II. FORMULATION 

In this Section we briefly discuss the basic equations, 
as well as their numerical implementation. For a more 
detailed discussion, we refer the reader to Section III of 
BSS for the gravitational fields (but point out some minor 
differences below) and to Section II of [34j for the hydro- 
statics. For corotating sequences using the Kerr-Schild 
metric, our notation will appear to be slightly different 
than that found in BSS, but the two sets of equations are 
completely equivalent, merely expressed in a different set 
of variables. 



A. Gravitational field equations 

The line element in 3+1 form is written as 

ds 2 — g^dx^dx 1 " 

= -a 2 dt 2 + 7 y (dx l + I3 i dt){dx j + I3' ] dt), (1) 

where a is the lapse function, (3 l the shift vector, 7^ the 
spatial metric, and g^ the space-time metric. The Ein- 
stein equations then split into two constraint equations 
- the Hamiltonian and the momentum constraint - and 
two evolution equations - one for the spatial metric 7ij 
and one for the extrinsic curvature 

Ka =~ (dtlij + Difr + Djfy) , (2) 

where Di denotes the covariant derivative associated with 
7ij . Using the conformal decomposition 

lij=i>%j, (3) 

where ip is the conformal factor and 7^ the background 
metric, the Hamiltonian constraint coupled with the 
trace part of the evolution equation of becomes 

D 2 a = Att^S + ~ip- 8 A i:j A ij + ^R+ \^ k2 
+t-0 l D i K - [{b lV )(D jV ) + (p i0 )(D ja )] ,(4) 
where we have defined 

cr = ln(m/; 2 ); (5) 
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this quantity was denoted u 0" in [Hj]. Di, D 2 — 
j^DiDj, Rij, and R are, respectively, the covariant 
derivative, covariant Laplace operator, Ricci tensor, and 
scalar curvature with respect to 7^ . We have also decom- 
posed the extrinsic curvature into its trace and traceless 
parts, 

K ij = ij- w A i3 + -7" if. (6) 

In the derivation of ([4]) we have assumed dt K — 0, which 
results in 

D 2 v = 4tt^ 4 (p + S)+ ^AijAv + \ip A K 2 

+^/3 i D i K-y j (D i v)(D j a) (7) 
a 

for the quantity 

v = lna. (8) 

Finally, we assume that dtfij = for quasiequilibrium 
configuration, so that ((2|) yields 

A» = ^ + &[? - p ij D k /3 k ^ . (9) 

Inserting this relation into the momentum constraint 
yields 

D 2 (3 l + \b l Djfr = 16ira^f - R)(3 j 
3 

4 ~ 2rv - ~ 

+ -aD l K- -pAVDj&r-k/). (10) 

The matter quantities on the right-hand side of equations 
(jl]), |Q and (fTC)|) are projections of the stress-energy ten- 
sor 

= (po + pi + P)u l± u u + Pg^, (11) 

where is the fluid 4-velocity, po the baryon rest-mass 
density, pi the internal energy density, and P the pres- 
sure. We then define 

p = n^T^, (12) 

f = -%n u T^, (13) 

S y = w^T"", (14) 

5 = 7 y '%, (15) 

where — (—a, 0,0,0) is the future-directed unit nor- 
mal vector. 



B. Background black hole metric 

We account for the neutron star's black hole compan- 
ion by choosing a background solution that represents a 
Schwarzschild metric. In BSS the Schwarzschild solution 



was expressed in Kerr-Schild coordinates. In this paper 
we analyze the effect of this choice by comparing with 
the Schwarzschild background expressed in isotropic co- 
ordinates. Since the two coordinate system represent two 
distinct slicings of the Schwarzschild metric, we solve the 
constraint equations on different spatial slices, and there- 
fore have no reason to expect that the resulting solutions 
to the constraint equations are physically identical. 

TABLE I: Lapse function «bh, shift vector /3 BH , confor- 
mal factor i/>bh, and the conformally related spatial metric 
jij for the Schwarzschild metric in Kerr-Schild (K-S) and 
isotropic (ISO) coordinates. Mbh is the black hole mass, 
ran = V X 2 + Y 2 + Z 2 the distance from the black hole cen- 
ter, and we define Hbh = A/bh/Vbh and U = f = X 1 /Vbh. 
Note that tbh in Kerr-Schild coordinates indicates a different 
displacement than one measured in isotropic ones. 





K-S 




ISO 


«BH 


(i+2 J ff BH r 1/2 


1 
1 


- Hbh/2 
+ Hbh/2 


Pbh 


2a| H _ffBH^ 1 







ipBH 


1 


1 


+ Hbh/2 


% 


Tjij + 2HBHklj 




Vij 



In Table U we list the background metric quantities for 
both Kerr-Schild (KS) and isotropic (ISO) coordinates 
(43| . From these metric quantities a number of back- 
ground quantities - for example i?* and K - that enter 
equations ([!]), ([7]) and (jTDJ) are derived. For Kerr-Schild 
coordinates these quantities can be found in BSS; for 
isotropic coordinates they are either zero or trivial, since 
the background metric "/ij is flat and the shift vanishes. 

To aid in the numerical solution we decompose the 
metric quantities into contributions from the neutron star 
and the black hole as follows. The total lapse function 
and conformal factor are decomposed into the product of 
a neutron star part and a black hole part (not a sum as 
in BSS), such that 

a = hnsobh, (16) 

1p = ^NS^BH- (17) 

The neutron star part is calculated by solving Poisson- 
like equations while the black hole part is given by the 
background metric, shown in Table HI The product form 
is required in order to decompose v and tr, defined by 
Eqs. ((HI) and ([5]), as a sum of the neutron star and black 
hole parts, such that 

V = I^NS + ^BH, (18) 

a = ons + obh, (19) 

where j/ N s = hra N s, ^bh = lnctBH, cns = MonsV'ns)' 
and <tbh = ln(aBH^BH)- The decomposition (|16p and 
(fT?)) are different from those employed by BSS, but for- 
mally equivalent. 
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For the shift vector, we decompose as 

^'=#IS+/%H + /&t> (20) 

where /3^ s and /3g H are the neutron star and black hole 
contributions to the shift vectors seen by the inertial ob- 
server, and /3* ot is the rotating shift vector defined as 

(3l ot = e^ k n,X k = Q(-Y, X, 0), (21) 

where £lj is the orbital angular velocity vector. We define 
the Z-axis to be parallel to the rotation axis, so that 
51* = (0, 0, f2) in Cartesian coordinates. 

As in BSS, we assume extreme mass ratios Mbh 3> 
A/nSj which simplifies the problem in several ways. In 
this limit, we may assume that the neutron star affects 
the space-time only in a region around the neutron star 
itself, so that we can restrict the computational domain 
to a neighborhood of the neutron star. This means that 
we do not need excise a black hole singularity from the 
numerical grid. We may also assume that the rotation 
axis coincides with the center of the black hole, which 
eliminates the need for an iteration to locate the axis of 
rotation (see Fig. [T]). 

For the isotropic background, the solution is symmet- 
ric both across the equatorial plane (i.e. the X—Y plane) 
and the X—Z plane. As discussed in BSS and [4J|, the 
presence of a non-vanishing trace of the extrinsic curva- 
ture eliminates the symmetry across the X—Z plane for 
a Kerr-Schild background, so that in the latter case we 
can assume a symmetry only across the equatorial plane. 

Evidently, the numerical implementation of the 
isotropic background is much easier than that of the 
Kerr-Schild background. However, isotropic coordinates 
have the disadvantage that the lapse aBH vanishes on the 
black hole horizon, which would cause problems in ([9]) 
when the black hole is included in the computational do- 
main. The same property also causes problems in dynam- 
ical simulations, when these quas ieq uilibrium models are 
used as initial data (thus, in [14j, Poisson-like equations 
are solved for i/j and aift, as in BSS). Kerr-Schild coordi- 
nates have the advantage that the coordinates smoothly 
extend into the black hole interior, which eliminates both 
of these problems. 

The above decompositions are then inserted into equa- 
tions dU) , and ([TO]) , which are then solved for the neu- 
tron star contributions. We list the resulting equations 
in Appendix [Al 

C. Hydrostatic equations 

For the cases considered here, the equations of rela- 
tivistic hydrodynamics reduce to the the integrated Eu- 
ler equation as well as the equation of continuity (which 
is satisfied identically for corotating stars). There are 
many sets of notation used to describe this formulation 
ES EE SI EE El, but we follow that of 0, with one 
important exception: we denote the adiabatic index of 




FIG. 1: Coordinate systems. The origin of the "global" co- 
ordinates (X, Y, Z) is located at the center of the black hole, 
and that of the "local" coordinates (x,y,z) at the point of 
maximum density within the neutron star. Here, we assume 
an extreme mass ratio Mbh Mns, so that the rotation 
axis of the orbit coincides with the center of the black hole. 
We define the Z-axis to point in the direction of the rotation 
axis. The relation between the global and local coordinates 
is given by X = x + Xns, Y = y, and Z — z, where Xns 
denotes the orbital separation between the black hole and the 
neutron star. 

our polytropic EOS by T and the Lorentz factor of the 
matter by 7, whereas their notation defines them the 
other way around. 

1. Integrated Euler equation 

The first integral of the Euler equation, common to 
both corotating and irrotational configurations, can be 
written as 

ha — = constant, (22) 
7o 

where h = (p + pi + P)/p is the fluid specific enthalpy. 
Here 7 is the Lorentz factor between the fluid and the 
co-orbiting observer and 70 that between the co-orbiting 
observer and the inertial frame. Here we define the "in- 
ertial" frame as the frame corresponding to normal ob- 
servers for whom the shift vector goes to zero at large 
tbh- For the co-orbiting frame, normal observers satisfy 
[3 l — > (ft x tbh) 1 at large rpH- A detailed derivation can 
be found, for example, in [341 ]. The Lorentz factors can 
be written as 

7 = 7n7o(l-7yC/ 1 ^), (23) 
70 = (l- 7 «^)- 1/2 , (24) 
7n = (l-^WW)- 1 / 2 , (25) 

where Uq is the orbital 3-velocity with respect to the 
inertial observer, given by Uq — f3 l /a, where the shift 
vector j3 l is measured by a co-orbiting observer. The 
quantity "f n denotes the Lorentz factor between the fluid 
and the inertial observer, and can also be expressed as 
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7 n = otu where u is the time component of the fluid 4- 
velocity u M . The quantity U % is the fluid 3- velocity with 
respect to the inertial observer (see Eq. (27) of |34|); for 
corotating binary systems, the fluid 3- velocity seen by the 
co-orbiting observer vanishes, we obtain U l = Uq. For 
irrotational binary systems, U l can expressed in terms of 
a velocity potential iff as 



(26) 



where Di is the covariant derivative with respect to 7^ . 
In this case, 7 n is written as 



7n = 1 



(27) 



Taking the logarithm of Eq. (|22p , we obtain the final 
form of the integrated Euler equation, 



-Sent + v — ln7o + In 7 = constant, 



(28) 



where -ff C nt = In h. For corotating binary systems we 
have 7=1 and hence In 7 = 0. 



2. Equation of continuity 

Having taken into account the helical symmetry, we 
rewrite the equation of continuity 



Tl / Tl 

-V^V A1 * + (V^ l *)V A1 (- 



0, 



(29) 



where n is the fluid baryon number density and V M the 
covariant derivative associated with t^,„, in the 3+1 form 

nD'DiV + {D l n){D^) = h^U^D.n 



{D l ^)D t In - + hW Q Dan + nhK ln . (30) 



Introducing an auxiliary quantity which represents the 
logarithmic derivative of enthalpy with respect to baryon 
number density, 



c 



dlniL 



dlnn 



(31) 



we can rewrite Eq. (|30|) as 



(H ent DH + fWjH^D^ = ip 4 h ln W DiH ont 



+CH C 



fW^D^H^t -a) + ^hU l Q b iln 



(32) 



We will show the final form of equations for the cases of 
the Kerr-Schild and isotropic coordinates in Appendices 
lATbl and [A2bl 



D. Equation of state 

We adopt a polytropic EOS for the neutron star of the 
form 



P = K P T a 



(33) 



where T denotes the adiabatic index and k is a constant. 
All results shown in this paper are for Y = 2. Since 
dimensions enter the problem only through k, it is con- 
venient to rescale all dimensional quantities with respect 
to the length scale 



-^poly — K ^ ^ \ 



(34) 



For later comparison with the results of BSS, we also 
introduce the dimensionless quantity 



P 

9= —, 
Po 



(35) 



in terms of which we can express the baryon rest-mass 
density and the pressure as 



p-2 V(r-i) 



P - R- 2 a r/{r - 1) 



E. Boundary condition 



(36) 
(37) 



For any Poisson-like elliptic equation of the form 

A<f> = s, (38) 

where s is a source term with compact support, we can 
express the exterior solution as 

00 / 

$(r,0,¥O=£ simY lm (8,<f)r-V +1 \ (39) 

i=0 m=— I 

Here, si m are the multipolc moments, defined as 



si, 



(40) 



and where Rb is the extent of the non-zero domain of the 
source term. The source terms of the Poisson-like equa- 
tions (fATj) - (JXHJ or (|A12j) - (|A14|) arc not compact in the 
situation we consider; however, since they fall off with a 
steep power-law dependence on the radius, ignoring the 
contribution from outside the computational domain in- 
troduces only minimal errors. In our code we match the 
numerical solution to (|39p at the outer boundaries, which 
we refer to as a multipole boundary condition. 

We truncate the expansion of the solution at a pre- 
determined value I = Z max , here setting Z max = 4 through- 
out. We set sim = for all terms with I > Z max . 

Multipole boundary conditions, which allow us to 
achieve much greater accuracy, are convenient in spectral 
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applications, but less so in finite difference implementa- 
tions. BSS, for example, used Robin boundary conditions 
which enforce that fields fall of with a certain power of 
1/r. In Sec. IIII A 21 we will compare the results of these 
two boundary conditions (and we also refer to [14j for a 
more detailed discussion of the boundary conditions). 



F. Determination of the orbital angular velocity 

To determine the orbital angular velocity, we require 
a force balance along the X-axis at the center of the 
neutron star, 



dX 



(x NS ,o,o) 



= 0, 



(41) 



where Xns is the AT-coordinate of the center of the 
neutron star relative to the black hole. Equation (|4"Tj) 
means that the neutron star has its maximum enthalpy 
(and thus density) at the position (X^g, 0, 0), which cor- 
responds to the origin of the local coordinate system 
[x, y, z) in which we solve the equations (see Fig. [1]). 
Inserting Eq. (|28]) into the condition (|4"Tj) . we obtain 



9 1 
~dX 



= ^ — \v + in 7 

(x N3 ,o,o) dX y " 



(X NS ,0,0) 



(42) 



Equation (|42|) is solved algebraically for the orbital angu- 
lar velocity fi, which enters through the Lorentz factors 
and the shift decomposition (|2"0|) on both sides of the 
equation. We leave the dependence implicit on the right- 
hand side, but write out the dependence explicitly on the 
left-hand side (see Appendices I A 1 cl and I A 2"c|) . 



Assuming extreme mass ratios allows us to neglect 
tidal effects of the neutron star on the black hole and 
to restrict the computational domain to a neighborhood 
of the neutron star. This means that we cannot evaluate 
the black hole's irreducible mass Afj rr . Instead, we keep 
the background mass Mbh constant in this paper. The 
difference between M- m and Mbh is of the order of the 
binary's binding energy (see footnote [35] in BSS), which 
is much smaller than either M; rr or Mbh in the limit of 
extreme mass ratios. 

Other global integrals that we might be interested in 
are the ADM mass and the angular momentum, and es- 
pecially their change along equilibrium sequences as a 
function of the binary separation. However, for the same 
reasons as explained in the previous paragraph, we can- 
not capture contributions of the binding energy to these 
quantities, meaning that we would find errors as large 
as the quantities that we are interested in. We there- 
fore postpone evaluation of these integrals until we have 
relaxed the assumption of extreme mass ratios. 



(d, 0, 0) 




(0, 0, 0) 



FIG. 2: Relation among the various neutron star "radii" r e , r a 
and Vb, as well as between the "orbital separations" Xns and 
d. r e is the half-diameter, which is defined as r e = (r a + rt) /2. 



G. Global Integrals 



It is reasonable to assume that during the binary inspi- 
ral both the neutron star's rest-mass and the black hole's 
irreducible mass are conserved. The rest mass is defined 



M, 



NS 



p u 



-gd 'fx, 



(43) 
(44) 



where g is the determinant of g^ u and 7 that of 7^ . The 
determinant 7 takes the form 

{1 + 2Hbh for the K— S background, 
1 for the flat background. 

In polytropic units (see Eq. (|34p ), we can normalize the 
baryon rest-mass in dimensionless form 



M NS 



M N s 



poly 



(46) 



H. Solution procedure 

We construct quasiequilibrium BHNS binaries in an 
iterative algorithm as follows. We start by preparing an 
initial guess, and then construct self-consistent solutions 
to the equations listed in Appendices IA II and IA 21 in an 
iteration that is similar to that of [34{. Since there are 
some differences in various steps of the process, however, 
it may be useful to describe the iteration in detail: 

1. Preparation of initial data for the main iteration 

(a) Construct a spherical star in equilibrium. 

(b) Set the orbital separation between the center 
of the neutron star and the black hole to X-^s ■ 
Here, the center of the neutron star is located 
at the maximum of the enthalpy (and thus the 
density). As we assume an extreme mass ra- 
tio, the absolute coordinate of the center of 
the black hole is (0,0,0), that of the neutron 
star (Xns, 0, 0); the rotation axis points in the 
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Z-direction, and goes through the origin of the 
absolute coordinate (see Fig.[T|). At this stage, 
we set the black hole mass to a given constant 
A/bh and never change it during the main it- 
eration that follows. 

2. Main computation 

(a) Introduce the black hole gravitational field as 
a fixed background. 

(b) Set the initial guess of the orbital angular ve- 
locity to the Keplerian frequency. 

(c) Solve the set of equations 

i. Gravitational field equations: Eqs. (|A1[) 

(|A3|) for Kerr-Schild background; 
Eqs. (|A12|) - (IA14j) for isotropic back- 
ground 

ii. Equation of continuity: Eq. (|A10[) for 
Kerr-Schild background; Eq. (|A15[) for 
isotropic background 

(d) Determine the new enthalpy field from the 
integrated Euler equation, Eq. ([28]) . and the 
new baryon rest-mass density from the neu- 
tron star EOS. 

(e) Calculate the orbital angular velocity from 
Eq. (|42[) by solving it algebraically. 

(f) Re-scale the radius of the neutron star and fit 
the outer boundary of the innermost domain 
to the surface of the neutron star. 

(g) Change the orbital separation to have the 
same ratio of X-^g/r e as the initially given 
value. Here r e is the half-diameter of the neu- 
tron star on the X-axis (see Fig. [5]). 

(h) Change the central enthalpy to fix the baryon 
rest-mass of the neutron star at a given value 
M NS . 

(i) Compare the new enthalpy field with that of 
old one, and check whether the relative dif- 
ference is smaller than the threshold or not. 

(j) If the condition ([2i]) is not satisfied, go back 
to (l2cl) and continue. 



The lack of a symmetry across the X— Z plane for 
a Kerr-Schild background introduces one more subtlety. 
The procedure we introduced above determines all the 
eigenvalues which appear in the present case, i.e., the or- 
bital angular velocity, the radius of the neutron star, and 
the integration constant of the Euler equation. However, 
the above procedure does not include a method to fix the 
position of the neutron star in the local coordinate system 
{x,y,z), in which we solve the field equations. Thanks 
to the X—Y plane symmetry (equatorial symmetry), the 
center of the neutron star is automatically located in the 
X—Y plane. Additionally, the procedure for determi- 
nation of the orbital angular velocity (f2"e)) fixes the local 



maximum of the neutron star rest-mass density along the 
X-axis at X = Ans (x = 0). Lacking an additional con- 
straint, the neutron star center can fall at any position 
in the ^-direction, so long as the A-coordinate takes the 
proper value. To define our configurations unambigu- 
ously, we must require that the position we impose in 
the procedure (|2ep is the global maximum of the neutron 
star rest-mass density. This requires the ^-derivative of 
the enthalpy to be zero as well as the A-derivative at the 
point (Ans, 0, 0). To do so, we first introduce a function 
f(y) = 1 — Ay, where A is a constant defined by 



A = 



1 dH e 



H, 



dy 



(X NS ,0,0) 



(47) 



When we multiply this function by the enthalpy and 
define a modified enthalpy -ff mo d = fH ent , the mod- 
ified enthalpy has its global maximum at the position 
(Ans, 0,0). During the iteration, this modified enthalpy 
term drags the neutron star to the proper position in 
the Y-direction. When the enthalpy maximum is prop- 
erly located on the A-axis, we recover A = and hence 
/ = 1. We insert this procedure between J2f} and ( |2"g| ) 
for the case of the Kerr-Schild background metric. 

On the other hand, we have a X—Z symmetry plane 
when using the isotropic background. Thus, the Y- 
derivative of the enthalpy always becomes zero on the 
A-axis, and the global maximum of the enthalpy must 
fall on the A-axis. 



III. CODE TESTS 

Our numerical code is based on the spectral methods 
libraries developed by the Meudon relativity group [35[, 
and have been tested, used, and found to be highly ac- 
curate in several previous app lications involving binary 
neutron stars 0, H H3, S H El • In the following 
we present some tests that verify our code for BHNS bi- 
naries. We present self-consistency tests, in which we 
explore the sensitivity of the results to changes in the 
position of the outer boundary, the type of boundary 
condition, and the number of collocation points used in 
the spectral method, and compare with analytical and 
previous numerical results. 



A. Self-consistency tests 

We present here the results of the self-consistency 
checks for a neutron star mass A/ns = 0.05. Our con- 
vergence criteria are chosen so that we iterate our solver 
until the relative difference in mass between the given 
value and our numerical result is less than one part in 
10~ 7 . In the following tests, we will present values of the 
results with five digits, although this overstates their true 
accuracy, in order to show the magnitude of the changes 
in these quantities. 
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1. Position of the outer boundary 

As long as the outer boundary of our computational 
domain is sufficiently far away from both the neutron 
star and the black hole, any physical quantities should 
be unaffected by the exact location of the outer bound- 
ary. Note that the outer boundary is located between 
the neutron star and the black hole in the present work 
(see Fig. [T|). To test this, we set the orbital separation 
to XNs/ r e = 10, and change the position of the outer 
boundary from R^/r = 4 to 9, where r denotes the ra- 
dius of a spherical neutron star with the same baryon 
rest-mass. In Table HH we show the orbital angular 
velocity QMbb., the half-diameter of the neutron star 
f e = r e /i? P oiy, and the maximum of the density quantity 
<3max- One can see from Table \H\ that the radius of the 
neutron star and the maximum of the density quantity 
converge to better than order 10 -5 for 5 < i?e/^o < 8, 
while the convergence of the orbital angular velocity is 
only of order 10 -3 for a Kerr-Schild background and 10 -4 
for an isotropic background. 

For a separation ATnsAe = 10 and Mns = 0.05, no 
position of the outer boundary larger than R^/r^ > 9.7 
is permitted because the boundary will overlap the black 
hole singularity. Even for a smaller value, the source 
terms of Eqs. JM) - CM]) or (fAT2|) - (fAT4|) do not de- 
crease sufficiently quickly with radius because of the con- 
tributions from the black hole background metric. Thus, 
we have to choose a position of the outer boundary for 
each orbital separation whose position is neither close to 
the neutron star nor to the black hole. 

The relative error of the orbital angular velocity listed 
in Table HT1 is larger for the Kerr-Schild background than 
for the isotropic background. This tendency results in 
part from the absence of the X—Z plane symmetry in 
the Kerr-Schild background. We also see that the or- 
bital angular velocity determined by using a larger outer 
boundary radius is slightly larger than that found using 
a smaller one. Because of this ambiguity, we will show 
only three significant digits in the final results. 



2. Type of boundary condition 

Next, we show in Table IIIII the results of a compari- 
son between our two different boundary conditions. The 
first is the multipolc condition described in Sec. Ill Ei and 
second is the condition in which only the leading-order, 
fall-off term is taken for each metric component. It is 
found that the relative difference in the orbital angular 
velocity is smaller than order 10 -3 for the Kerr-Schild 
background and 10 -4 for the isotropic background. The 
leading-order fall-off boundary condition seems to overes- 
timate the angular velocity relative to the multipole con- 
dition. However, the relative error introduced by using 
different boundary conditions is an order of magnitude 
less than that introduced by changing the position of the 
outer boundary. 



3. Number of collocation points 

In the last self-consistency check we verify that phys- 
ical quantities are largely independent of the number of 
collocation points N r x Ng x N v used in spectral field 
solver. Here, N r , Ng, and N v denote the number of col- 
location points in the radial, polar, and azimuthal direc- 
tions, respectively. The orbital angular velocity seems to 
converge at a level of 10~ 4 in the Kerr-Schild background 
and < 10~ 5 in the isotropic background, for any number 
of collocation points larger than 25 x 17 x 16. Since the 
error in the orbital angular velocity is smaller than that 
induced by changing the position of the outer boundary 
for this number of points, we use a fixed number of col- 
location points N r x Ng x N v = 25 x 17 x 16 to save 
computational time. 



B. Comparison with previous results 

To date, previous attempts to model equilibrium 
BHNS binary sequences have included the construction 
of a corotating relativistic star in a black hole background 
metric [31], a corotating Newtonian star in a black hole 
background metric computed in the tidal approximation 
[HI, and a corotating or irrotational BHNS binary sys- 
tem in Newtonian gravity p9j . We now compare our 
relativistic results with these calculations. 



1. Comparison with BSS 

BSS provide results for a constant rest-mass sequence 
of corotating BHNS binaries in a Kerr-Schild back- 
ground, with Mns = 0.05 and Mbh/Mns = 10. We com- 
pare our results at two of the orbital separations listed in 
their Table I, d/r e — 8 and 5. Here d is the coordinate 
separation from the center of the black hole to the half- 
diameter of the neutron star on the X-axis, a quantity 
similar to their "— xbh" ■ We define 

d = X NS + ^A (48) 

where r a is the radius of the neutron star in the direction 
toward the black hole and that to the opposite side (see 
Fig. [J) . Since our definition of the orbital separation Xns 
is different from that found in BSS, we have to transform 
our separation Xns to d. However, since the relative 
difference between these two separations is less than 1% 
even for d/r e = 5, our results are insensitive to the change 
up to three significant digits. 

In Table [V] we present the results of the comparison. 
Note that we show the results at d/r e = 5 in the table 
here, but our results at this separation do not converge as 
well as our other results, only to the level of several parts 
in 10 -5 in the relative error of the enthalpy between two 
successive iteration steps. Thus we do not include those 
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TABLE II: Convergence test for the position of the outer boundary Rb ■ 



Kerr-Schild: M NS = 0.05, Xns/r e = 10.0, f = 1.1289 



Rs/ro 



Irrotation 

ymax 



Corotation 

UMbH fe Umax 



9.7328(-3) 1.0964 
9.7432(-3) 1.0962 
9.7515(-3) 1.0962 
9.7603(-3) 1.0962 
9.7760(-3) 1.0962 
9.9358(-3) 1.0964 



0.023493 
0.023493 
0.023493 
0.023493 
0.023493 
0.023488 



9.6399(-3) 1.1030 0.023339 

9.6500(-3) 1.1029 0.023338 

9.6580(-3) 1.1029 0.023338 

9.6664(-3) 1.1029 0.023338 

9.6807(-3) 1.1029 0.023337 

9.7890(-3) 1.1032 0.023330 



Isotropic: Mns = 0-05, X N s/r e = 10.0, fg = 1.1289 



Rs/ro 



Irrotation 



*?max 



Corotation 



?max 



9.1059(-3) 1.0960 0.023507 

9.1093(-3) 1.0960 0.023507 

9.1114(-3) 1.0960 0.023507 

9.1129(-3) 1.0960 0.023507 

9.1139(-3) 1.0960 0.023507 

9.1149(-3) 1.0960 0.023507 



9.0364(-3) 1.1017 0.023376 

9.0397(-3) 1.1017 0.023376 

9.0418(-3) 1.1017 0.023376 

9.0432(-3) 1.1017 0.023376 

9.0443(-3) 1.1017 0.023376 

9.0452(-3) 1.1017 0.023376 



TABLE III: Convergence test for the type of boundary condition. 



Kerr-Schild: M N 



0.05 









Irrotation 


Corotation 


Boundary condition Xns/fe 


-Rb /ro 


HMbh r e g max 


OMbh 


Qmax 


multipole 


10.0 


8 


9.7760(-3) 1.0962 0.023493 


9.6807(-3) 


1.1029 0.023337 


leading order 


10.0 


8 


9.7799(-3) 1.0962 0.023493 


9.6843(-3) 


1.1029 0.023337 


multipole 


8.0 


7 


1.3588(-2) 1.1062 0.023429 


1.3325(-2) 


1.1186 0.023126 


leading order 


8.0 


7 


1.3627(-2) 1.1064 0.023427 


1.3350(-2) 


1.1188 0.023125 


multipole 


6.0 


5 


1.9457(-2) 1.1583 0.023193 


1.8837(-2) 


1.1818 0.022558 


leading order 


6.0 


5 


1.9478(-2) 1.1584 0.023192 


1.8855(-2) 


1.1817 0.022558 


Isotropic: Mns = 0.05 








Irrotation 


Corotation 


Boundary condition ^ns/^b 


-Rb /ro 


QMbh f e (J max 


OMbh 


fe Qm ax 


multipole 


10.0 


8 


9.1139(-3) 1.0960 0.023507 


9.0443(-3) 


1.1017 0.023376 


leading order 


10.0 


8 


9.1152(-3) 1.0960 0.023507 


9.0455(-3) 


1.1017 0.023376 


multipole 


8.0 


7 


1.2430(-2) 1.1021 0.023467 


1.2261(-2) 


1.1126 0.023221 


leading order 


8.0 


7 


1.2432 (-2) 1.1022 0.023467 


1.2263(-2) 


1.1126 0.023221 


multipole 


6.0 


5 


1.7811(-2) 1.1381 0.023315 


1.7372(-2) 


1.1580 0.022800 


leading order 


6.0 


5 


1.7816(-2) 1.1382 0.023315 


1.7377(-2) 


1.1581 0.022800 



results in Sec. IIV A[ but show them here as a code test 
only. 



Our results agree with those of BSS within 2%, even 
though we calculate them using a completely different 
numerical method. 



2. Comparison with the results by Ishii et al. 129] 



In [28(, the critical value of a quantity 



VMrh/V r n J ' 



(49) 
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TABLE IV: Convergence test for the number of collocation points used in the spectral method field solver, N r x Ng x N v . 



Kerr-Schild: M NS = 0.05, X NS /r e = 10.0, R B /r = 8 



N r xNgX N v 


Irrotation 
OA/bh f e g max 


Corotation 
OMbh f e g max 


13 x 9 x 8 
17 x 13 x 12 
25 x 17 x 16 
33 x 21 x 20 
37 x 25 x 24 


9.8004(-3) 1.0961 0.023495 
9.7807(-3) 1.0962 0.023493 
9.7760(-3) 1.0962 0.023493 
9.7742(-3) 1.0962 0.023493 
9.7735(-3) 1.0962 0.023493 


9.7019(-3) 1.1030 0.023338 
9.6849(-3) 1.1029 0.023337 
9.6807(-3) 1.1029 0.023337 
9.6792(-3) 1.1029 0.023337 
9.6785(-3) 1.1029 0.023337 


Isotropic: M NS = 0.05, X NS /r e = 10.0, R B /r = 8 


N r X No X jV<p 


Irrotation 

fiMBH f e (Jmax 


Corotation 
OMbh fe g max 


13 x 9 x 8 
17 x 13 x 12 
25 x 17 x 16 
33 x 21 x 20 
37 x 25 x 24 


9.1144(-3) 1.0959 0.023511 
9.1140(-3) 1.0960 0.023507 
9.1139(-3) 1.0960 0.023507 
9.1139(-3) 1.0960 0.023507 
9.1139(-3) 1.0960 0.023507 


9.0468(-3) 1.1016 0.023379 
9.0443(-3) 1.1017 0.023376 
9.0443(-3) 1.1017 0.023376 
9.0443(-3) 1.1017 0.023376 
9.0443(-3) 1.1017 0.023376 



TABLE V: Comparison with the results of BSS, which are 
calculated in the Kerr-Schild background metric for corotat- 
ing BHNS binaries. The neutron star has a baryon rest-mass 
Mns = 0.05, and the mass ratio is Mbh/A/ns = 10. 





d/re 


d/MBH 


JJMbh 


fe 


Qmax 


BSS 


8.0 


17.8 


0.0133 


1.11 


0.0235 


Our results 


8.0 


17.9 


0.0133 


1.12 


0.0231 


BSS 


5.0 


12.6 


0.0223 


1.26 


0.0223 


Our results 


5.0 


12.7 


0.0222 


1.28 


0.0219 



was calculated for corotating BHNS binaries in Fermi 
normal coordinates. Here ro is the radius that a neutron 
star with the same mass would have in isolation. In these 
calculations, the neutron star is treated as a Newtonian 
star in the black hole background which is expanded up 
to fourth order in the parameter ro/A^s- The quantity 
above is defined such that if /i < fi C nt, the neutron star 
will be tidally disrupted by the black hole. By fitting 
their numerical results, they give an approximate formula 
for /i cr it which holds for separations ATns > 6Mbh with 
good accuracy. The formula, Eq. (186) of [28| . is given 
by 

Merits 14.9(l + 0.80^-) (50) 

for polytropic neutron star models with an adiabatic in- 
dex T = 2. 

We show the results of the comparison in Fig. [3] 
The thick solid and dashed lines correspond to sequences 
of corotating BHNS binaries with neutron star masses 
M~ns = 0.05 and 0.1, respectively, while the filled circle 



and square show the critical values for these sequences. 
These points are calculated by inserting the closest sep- 
aration we obtained for the corotating state in the Kerr- 
Schild background case into Eq. ([4"9"|) . The thin dotted 
and dot-dashed lines denote the critical value given by 
Eq. ([50]) for neutron star masses Mns = 0.05 and 0.1, 
where we use the radius of a spherical neutron star with 
the same baryon rest-mass; ro = 1.1289 for Mns = 0.05 
and fo = 0.98972 for Mns =0.1. for drawing all critical 
points. 

Evidently, our sequences do not quite reach the limits 
found in [28|]. This may be because our spectral method 
is unable to treat the cusp-like features that appear at the 
inner Lagrange point when the neutron star fills out its 
Roche lobe. We can therefore approach the tidal break- 
up separation only up to a few percent, which may ex- 
plain why our results for ^ cr ;t are somewhat larger than 
those found in [28]. 



3. Comparison with the results by Uryu & Eriguchi J2W 

In this section, we compare our results with those 
of as was done by BSS with good agreement. In 
[29I ] sequences of BHNS binaries in Newtonian gravity 
were computed for both irrotational flow and corota- 
tion. For comparison, we select data with A^ = 1.0 and 
A/s/A/bh = 0.1 in their Table 2 for the corotating case, 
and in Table 4 for the irrotational case. Here, A" denotes 
the polytropic index (N = l/QT — 1)) and Ms is the mass 
of the star (corresponding to Mns here). We present the 
results of a comparison of the orbital angular velocity 
in Table IVII There, we compare our data for a neutron 
star mass Mns — 0.05 and mass ratio Mbh/Mns = 10 
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Comparison with Ishii et al. 



64 
56 
48 
40 
"32- 
24- 
16 



~r 



— Sequence for M NS =0.05 

- ■ Sequence for M NS =0.1 

a , for M" =0.05 by Ishii et al. 

fcnt NS J 

— - 11 . for M. =0.1 by Ishii et al. 

^cnt NS J 



t Our 11 , forM.„=0.1 



crit NS 



'Oura . for M =0.05 

r rnt NS 



12 

X NS /M BH 



16 



20 



FIG. 3: Critical value of the tidal disruption parameter /j as a 
function of the orbital separation. The thick solid and dashed 
lines represent sequences of corotating BHNS binaries with 
neutron star masses Mns = 0.05 and 0.1, respectively; the 
filled circle and square show the critical values for both. The 
thin dotted and dot-dashed lines denote the critical values 
given by Eq. (|50p for neutron star masses Mns = 0.05 and 
0.1. 



TABLE VI: Comparison with the results of Uryu and Eriguchi 
[2!|. The results are compared for a neutron star mass Mns = 
0.05 with mass ratio Mbh/Mms = 10. Here, we use the radius 
of a spherical star, fo = 1.1289, to convert their results into 
our units. K-S and ISO denote the respective backgrounds. 





Irrotation 


Corotation 




flM BH 


fiM BH 


Uryu & Eriguchi 


0.0262 


0.0264 


Our results (K-S) 


0.0224 


0.0216 


Our results (ISO) 


0.0224 


0.0217 



in the Kerr-Schild and isotropic backgrounds. Since the 
sequences in [29| were computed using Newtonian grav- 
ity while we use full general relativity, we can compare 
only an invariant quantity, e.g., the orbital angular veloc- 
ity. In order to convert their Newtonian results to values 
comparable to our relativistic ones, we scale our results 
by the radius of a relativistic spherical neutron star with 
mass Mns = 0.05, fo = 1.1289. Our results are shown 
at the closest separation we calculate, which is several 
percent further away in distance from the real tidal dis- 
ruption point, because the spectral method is unable to 
treat cusp-like figures. Thus we find our calculated values 
are about 15 % smaller than those found in 29], but given 
that our configurations have separations several percent 
away from the real tidal disruption point, we can say that 
the agreement with the angular velocities in Table IVII is 



reasonable. 



C. Comparison with the relativistic Keplerian 
velocity in the limit of extreme mass ratio 

In the limit of an extreme mass ratio (Mbh/Mns 
1), the neutron star orbits on a circular, test-particle 
geodesi c. In thi s limit, we should recover Kepler's law, 
^ = a/A/eh/?" 3 , where r is the areal radius. In Kerr- 
Schild coordinates, the radial coordinate is areal, but in 
isotropic coordinates we have to transform according to 



(«) 2 nso = (l + 



M BH \ 2 



2ns< 



nso, 



(51) 



where nso denotes the isotropic radius. In Fig. [H we 
present the results for a mass ratio Mbh/A?ns = 100 with 
neutron star masses of A/ns = 0.05, 0.1, and 0.15. It is 
obvious that our results are in good agreement with the 
Keplerian frequency for both irrotational and corotating 
flows. 



IV. NUMERICAL RESULTS 

All results presented in this Section are for a polytropic 
index T = 2. The computational grid is divided into four 
domains (see Fig. [!}, each one of which is covered by 
N r x Ng x N v — 25 x 17 x 16 collocation points. We 



focus on a mass ratio of Mbh/Mns 
different sequences including 



10 and present 12 



• Kerr-Schild backgrounds and isotropic back- 
grounds 

• irrotational and corotational fluid flow 

• neutron star mass of Mns = 0.05, 0.1, and 0.15. 



The results are summarized in Tables IVTlI - IIXI for Kerr- 
Schild backgrounds and [X] - IXIII for isotropic back- 
grounds. Note that an isolated neutron star with the 
adopted T — 2 polytropic EOS has a maximum mass 

A/NS,max = 0.18. 



A. The Kerr-Schild background 

In Figs. [5] - [71 we show contours of the neutron star 
baryon rest-mass density profile at the closest separa- 
tion we calculate for neutron star masses Mns = 0.05, 
0.1, and 0.15, respectively. In each figure, the left panel 
shows the irrotational case and the right panel the coro- 
tating one. These figures correspond to the final lines in 
Tables I VIII - 11X1 Since the total shift vector induces an 
asymmetry in the neutron star with respect to the X—Z 
plane for a Kerr-Schild background, one can see a visible 
tilt in the stellar figure. 
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FIG. 4: The orbital angular velocity as a function of the orbital separation for a mass ratio Mbh/Mns = 100. The left and 
right panels show Kerr-Schild and isotropic backgrounds, respectively. In each panel, we present six different sequences: two 
choices for the rotational state (irrotation and corotation) for each of three different neutron star masses (Mns = 0.05, 0.1, 
and 0.15). The solid line in each panel denotes the frequency of a Keplerian orbit. Note that the horizontal axis shows the X 
coordinate in Kerr-Schild coordinates in the left panel, but isotropic coordinates on the right. 



TABLE VII: Physical parameters along a sequence with Mbh/Mns = 10 in Kerr-Schild backgrounds. The orbital separation 
Xns, orbital angular velocity il, half-diameter along the X-axis f e , maximum of the density parameter g m ax, relative change 
in the central energy density Se, mass-shedding indicator x> an d radius of the outer boundary Rb are shown. The baryon 
rest-mass of the neutron star is Mns = 0.05. M§ s denotes the gravitational mass of a spherical neutron star with the same 
baryon rest-mass. M§ s /ro is the compaction parameter of an isolated neutron star. 



Kerr-Schild: M NS = 0.05, M§ s /r = 0.0415, fg = 1.1289 
Irrotation 



Xns / r e 


Xns/Mbh 


CIMbh 


fe 


*?max 


5e 


X 


Rb/to 


10.0 


21.9 


9.78(-3) 


1.10 


0.0235 


-2.41(-3) 


0.977 


8 


8.0 


17.7 


1.36(-2) 


1.11 


0.0234 


-5.16(-3) 


0.908 


7 


7.0 


15.7 


1.62(-2) 


1.12 


0.0234 


-8.39(-3) 


0.834 


6 


6.0 


13.9 


1.95(-2) 


1.16 


0.0232 


-1.54(-2) 


0.699 


5 


5.2 


12.8 


2.24(-2) 


1.23 


0.0229 


-2.89(-2) 


0.505 


5 








Corotation 








X NS /r e 


X^s/Mbh 


fiM BH 


fe 


*?max 


5e 


X 


-Rb /ro 


10.0 


22.1 


9.68(-3) 


1.10 


0.0233 


-9.17(-3) 


0.957 


8 


8.0 


17.9 


1.33(-2) 


1.12 


0.0231 


-1.83(-2) 


0.875 


7 


7.0 


15.9 


1.58(-2) 


1.14 


0.0229 


-2.73(-2) 


0.793 


6 


6.0 


14.2 


1.88(-2) 


1.18 


0.0226 


-4.30(-2) 


0.654 


5 


5.2 


13.0 


2.16(-2) 


1.25 


0.0220 


-6.51(-2) 


0.460 


5 



In the left panel of Fig. [5j one can see a slight oscilla- 
tion on the stellar surface. During the computation, we 
fit the boundary of the innermost domain to the stellar 
surface, allowing us to solve equilibrium configurations 
accurately, particularly in irrotational cases. Once a cusp 
develops on the stellar surface, it becomes impossible to 
adapt the innermost domain to the stellar surface. Using 
spectral methods, we express all quantities by summation 



over a finite number of differentiable functions. If we ap- 
ply this method to such cusp- like figures, large numerical 
errors are induced. This situation is known as the Gibbs 
phenomenon. Since the stellar surface becomes highly 
distorted for very close configurations, even prior to the 
appearance of a cusp, we have to stop the sequence when 
such oscillations appear. 

In order to investigate how close our results come to the 
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TABLE VIII: Same as Table IVTT1 but for M NS = 0.1. 



Kerr-Schild: M NS = 0.1, M§ s /r = 0.0879, fg = 0.98972 









Irrotation 








s*. IN b / ' e 


X mo / M//du 
V IN b / JVJri±l 


li L J. H X3 J^l 


r e 


Qmax 


Se 


v 

A 


R / 

JIB / TQ 


14.0 


12.8 


2.19(-2) 


0.915 


0.0585 


-6.69(-3) 


1.06 


8 


12.0 


10.9 


2.82(-2) 


0.905 


0.0582 


-1.10(-2) 


1.05 


8 


10.0 


8.96 


3.78(-2) 


0.896 


0.0577 


-2.06(-2) 


1.02 


7 


8.0 


7.20 


5.29(-2) 


0.901 


0.0564 


-4.31(-2) 


0.912 


5.5 


7.4 


6.78 


5.88(-2) 


0.916 


0.0554 


-6.16(-2) 


0.821 


5.5 








Corotation 








^Ns/r e 


Ans/Mbh 


QMbu 


fe 


9max 


Se 


X 


Rb /ro 


14.0 


12.8 


2.18(-2) 


0.918 


0.0582 


-1.07(-2) 


1.05 


8 


12.0 


10.9 


2.79(-2) 


0.910 


0.0578 


-1.76(-2) 


1.03 


8 


10.0 


9.06 


3.71(-2) 


0.906 


0.0570 


-3.24(-2) 


0.990 


7 


8.0 


7.36 


5.15(-2) 


0.920 


0.0550 


-6.87(-2) 


0.860 


6 


6.8 


6.54 


6.13(-2) 


0.961 


0.0527 


-1.09(-1) 


0.684 


5 



TABLE IX: Same as Table jVUl but for M NS =0.15. 



Kerr-Schild: M NS = 0.15, M§ s /r = 0.145, f = 0.81526 
Irrotation 



X NS /r e 


Xns/Mbh 




fe 


*?max 


Se 


X 


-Rb /ro 


22.0 


10.7 


2.84(-2) 


0.733 


0.126 


-5.76(-3) 


1.10 


8 


20.0 


9.65 


3.34(-2) 


0.724 


0.126 


-8.32(-3) 


1.11 


8 


18.0 


8.55 


4.02(-2) 


0.712 


0.125 


-1.26(-2) 


1.13 


8 


16.0 


7.45 


4.97(-2) 


0.698 


0.124 


-2.09(-2) 


1.15 


8 


14.4 


6.57 


6.05(-2) 


0.683 


0.123 


-3.04(-2) 


1.16 


8 








Corotation 








-Xns / r e 


Ans/A^bh 


J2M BH 


fe 


*?max 


Se 


X 


Rb/to 


22.0 


10.8 


2.83(-2) 


0.734 


0.126 


-7.55(-3) 


1.10 


8 


20.0 


9.67 


3.32(-2) 


0.725 


0.125 


-1.09(-2) 


1.11 


8 


18.0 


8.58 


3.99(-2) 


0.715 


0.125 


-1.64(-2) 


1.12 


8 


16.0 


7.49 


4.91 (-2) 


0.702 


0.124 


-2.68(-2) 


1.13 


8 


14.0 


6.41 


6.27(-2) 


0.686 


0.121 


-4.86(-2) 


1.14 


8 



proper tidal disruption points, we introduce a sensitive 
numerical indicator for the mass-shedding point, defined 
as 

COmp /rrt\ 

X= (dH ont /dr) po i c ' 1 ' 

where the numerator denotes the radial derivative of en- 
thalpy of the neutron star toward the companion star in 
the equatorial plane on the X-axis and the denominator 
that toward the pole of the neutron star. Both terms 
are evaluated on the stellar surface. This indicator takes 
the value unity at infinite orbital separation, and goes to 
zero at the mass-shedding point. 



We show numerical results for \ in Fig. [8] The solid 
line with circles represents the irrotational case and the 
dashed line with squares the corotating one. As explained 
above, it is difficult for the spectral method to treat cusp- 
like figures, so we stop the calculation of sequences before 
reaching x = 0, which would correspond to the mass- 
shedding point. However, we can extrapolate our results 
in polynomial functions to predict the orbital separation 
at the tidal disruption point, \ = 0- The extrapolations 
are shown as a dotted line for the irrotational case and 
a dot-dashed line for the corotating case in Figs. [5Ja) 
and (b). For a neutron star of mass Mns = 0.15, we 
do not attempt to extrapolate our results to the tidal 
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Kerr— Schild (Irrotation) 



Kerr— Schild (Corotation) 
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FIG. 5: Baryon rest-mass density contours for a neutron star with Mns = 0.05 in the Kerr-Schild background, at the point of 
closest binary separation along our sequence. The left and right panels show the irrotational and corotating cases, respectively. 
The position of the density maximum is Xns/Mbh = 12.8 for the irrotational case and -Xns / Mbh = 13.0 for the corotating 
one. 
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FIG. 6: Same as Fig. [5j but for Mns = 0.1. The position of the density maximum is Xns/Mbh = 6.78 for the irrotational 
case and Xns/Mbh = 6.54 for the corotating one. 



disruption limit because the closest separations we can 
calculate are too far away, for both rotation states, to 
make a reasonable fit. Clearly the extrapolation intro- 
duces considerable error, so that we are able to predict 
the tidal separation only with modest accuracy. 

We find that the state of rotation has a small effect of a 
few percent on the tidal separation, which is in agreement 
with earlier findings (e.g. 27]). For Mns = 0-05 we 
find approximately ryd = 11.5Mbh (which is in good 
agreement with the value 11.9Mbh found in BSS), and 
for Mns = 0.1 we find approximately rtid = 5Mbh- 



Finally, we examine the relative change in the central 
energy density of the neutron star to that of an isolated, 
spherical neutron star with the same baryon rest-mass. 
Here, the total energy density and the relative change of 
its central value are, respectively, defined by e = po + pi 
and 

' ' (53) 



5e = 



1, 



e c ,o 



where e c is the central energy density of the neutron star 
and e C) o that of a spherical neutron star with the same 
baryon rest-mass. We see from Tables I VIII - IIXI that 
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Kerr— Schild (Irrotation) 



Kerr— Schild (Corotation) 



in 
o 



pq 



m 
o 




6.5 

X / M 



o 



pq 



d - 




BH 



6 6.5 

X / M BH 



FIG. 7: Same as Fig. [5] but for Mns = 0.15. The position of the density maximum is Xj-js/Mbh = 6.57 for the irrotational 
case and Xns/Mbh = 6.41 for the corotating one. 
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FIG. 8: Mass-shedding indicator \i defined by Eq. (|52p. as a function of the orbital separation for the neutron star mass 
(a)MNS = 0.05, (b)MNS = 0.1 and (c)Mns = 0.15, in Kerr-Schild backgrounds. Mass-shedding will begin when \ — 0- The 
thick solid line with circles represents the irrotational case and the thick dashed line with squares the corotating one. The thin 
dotted and dot-dashed lines are extrapolations for the results of the irrotational and corotating cases, respectively. 



the total energy density decreases as the orbital sepa- 
ration decreases for all sequences. At fixed orbital sep- 
aration and neutron star mass, the decrease in the to- 
tal energy density is larger for the corotating case than 
for the irrotational one. We attribute this effect to the 
neutron star spin, since the rotation expands the star's 
volume and results in a larger decrease in the central 
energy density. Roughly speaking, the amount of rela- 
tive decrease of the central energy density at the mass- 
shedding point is ~ 5% for an irrotational neutron star 
of mass Mns = 0.05, and ~ 10% for a corotating one, if 
we extrapolate our results. For a neutron star of mass 
Mns = 0.1, the relative decrease is about ~ 15% for 
the irrotational case and about ~ 20% for the corotat- 
ing one. This implies that a more massive neutron star 
undergoes a larger decrease in the energy density than a 
less massive one. We note that the central energy density 



decreases monotonically, even when the half-diameter on 
the X-axis, f e , decreases, as we find for sequences with 
Mns = 0.1 and 0.15. This occurs, in part, because f e 
is measured as a coordinate length, not as the proper 
length. Also, the diameter along the X-axis is not the 
primary axis in Kerr-Schild coordinates, because the stel- 
lar configuration is tilted due to asymmetries. 



B. The isotropic background 

In Figs. [3] - [Til we show contours of the baryon 
rest-mass density for neutron stars in an isotropic back- 
ground. As in the Kerr-Schild background, the left panel 
is for the irrotational case and the right panel for the 
corotating one. These figures are taken from the config- 
urations shown on the last lines of Tables |X] - IXIH the 
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TABLE X: Same as Table [VTTl but for isotropic backgrounds. 



Isotropic: M NS = 0.05, M§ s /r = 0.0415, f = 1.1289 
Irrotation 



ANS/r e 


V / 71 /T 

Ans/A^bh 




r e 


*?max 


oe 


X 


UB/ro 


10.0 


21.9 


9.11(-3) 


1.10 


0.0235 


-1.78(-3) 


0.942 


8 


8.0 


17.6 


1.24(-2) 


1.10 


0.0235 


-3.52(-3) 


0.883 


7 


6.0 


13.7 


1.78(-2) 


1.14 


0.0233 


-1.01(-2) 


0.719 


5 


5.0 


12.0 


2.13(-2) 


1.20 


0.0231 


-2.08(-2) 


0.514 


4.5 


4.7 


11.6 


2.24(-2) 


1.23 


0.0230 


-2.50(-2) 


0.406 


4 








Corotation 










Ans/Mbh 


OMbh 


r e 


Qm&x 


Se 


X 


Rs/ro 


10.0 


22.0 


9.04(-3) 


1.10 


0.0234 


-7.49(-3) 


0.925 


8 


8.0 


17.8 


1.23(-2) 


1.11 


0.0232 


-1.42(-2) 


0.855 


7 


6.0 


13.9 


1.74(-2) 


1.16 


0.0228 


-3.25(-2) 


0.681 


5 


5.0 


12.2 


2.07(-2) 


1.22 


0.0223 


-5.35(-2) 


0.495 


4.5 


4.7 


11.8 


2.17(-2) 


1.26 


0.0221 


-6.21(-2) 


0.413 


4 



TABLE XI: Same as Table H but for Af NS = 0.1 



Isotropic: M N s = 0-1, M§ s /r = 0.0879, f = 0.98972 









Irrotation 








A N s / r e 




OMbh 




?max 


Se 


X 


Rb/to 


14.0 


12.9 


1.94(-2) 


0.919 


0.0586 


-3.62(-3) 


0.980 


8 


10.0 


9.01 


3.15(-2) 


0.901 


0.0583 


-8.74(-3) 


0.936 


7 


8.0 


7.16 


4.27(-2) 


0.895 


0.0579 


-1.66(-2) 


0.863 


6 


6.0 


5.51 


5.96(-2) 


0.919 


0.0568 


-3.70(-2) 


0.659 


4.5 


5.6 


5.24 


6.35(-2) 


0.936 


0.0563 


-4.55(-2) 


0.577 


4.5 








Corotation 








A N s/r e 


Ans/A^bh 


nM B H 


fe 


Qraax 


Se 


X 


-Rb /ro 


14.0 


12.9 


1.93(-2) 


0.921 


0.0585 


-6.62(-3) 


0.973 


8 


10.0 


9.07 


3.12(-2) 


0.907 


0.0579 


-1.67(-2) 


0.918 


7 


8.0 


7.25 


4.20(-2) 


0.907 


0.0571 


-3.11(-2) 


0.837 


6 


6.0 


5.65 


5.78(-2) 


0.942 


0.0550 


-6.79(-2) 


0.634 


5 


5.6 


5.37 


6.16(-2) 


0.958 


0.0545 


-7.81(-2) 


0.564 


4.5 



closest separation we compute for each sequence. One 
can clearly see that there exists a symmetry with respect 
to the X—Z plane in addition to the X—Y plane symme- 
try (equatorial plane symmetry) in isotropic background. 
The half-diameter of the neutron star along the X-axis, 
f e , should be elongated by the tidal force from the black 
hole, but as in the Kerr-Schild background it decreases as 
the orbital separation decreases for Mns = 0.1 and 0.15, 
and then increases at closer separations for Mns = 0.1 
because of coordinate effects. 

The indicator for mass-shedding from the neutron star, 
X, is shown in Fig. [T2] Solid lines with circles and 
dashed lines with squares in each panel represent the 



irrotational case and corotating one, respectively. The 
dotted lines are extrapolations of our data for the irrota- 
tional case and dot-dashed lines those for the corotating 
case shown in Figs. [12] (a) and (b). For isotropic back- 
grounds, x does not exceed unity even for high compact- 
ness, and decreases monotonically when the orbital sep- 
aration decreases. From the extrapolations, we see that 
the orbital separations at the mass-shedding point are 
X_ NS ~ 10.5Af B H for M NS = 0.05 and AT N s ~ 4M BH for 
Mns = 0.1 with an error of approximately 10%. These 
values are coordinate separations, and should therefore 
not be compared immediately with the results in a Kerr- 
Schild background. We will compare our findings in much 
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TABLE XII: Same as Table fx] but for M NS = 0.15 



Isotropic: M N s = 0-15, M§ s /r = 0.145, f = 0.81526 
Irrotation 



X]s!s/re 


Xns/Mbh 


QMbh 


r E 


(Jmax 


Se 


X 


Ra/ro 


20.0 


9.75 


2.83(-2) 


0.732 


0.126 


-5.24(-3) 


0.995 


8 


16.0 


7.59 


3.96(-2) 


0.712 


0.126 


-9.67(-3) 


0.988 


8 


12.0 


5.46 


6.06(-2) 


0.682 


0.124 


-2.31(-2) 


0.959 


7 


10.0 


4.43 


7.82(-2) 


0.665 


0.122 


-4.22(-2) 


0.913 


6 


9.0 


3.96 


8.94(-2) 


0.659 


0.120 


-6.06(-2) 


0.868 


5.5 








Corotation 








X NS /r e 




fiM BH 


fe 


*?max 


Se 


X 


Rb /ro 


20.0 


9.77 


2.82(-2) 


0.733 


0.126 


-7.23(-3) 


0.992 


8 


16.0 


7.61 


3.94(-2) 


0.714 


0.125 


-1.35(-2) 


0.982 


8 


12.0 


5.50 


6.00(-2) 


0.687 


0.123 


-3.18(-2) 


0.948 


7 


10.0 


4.49 


7.68(-2) 


0.674 


0.120 


-5.60(-2) 


0.899 


6 


9.4 


4.20 


8.30(-2) 


0.671 


0.119 


-6.74(-2) 


0.874 


5.5 
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FIG. 9: Baryon rest-mass density contours for neutron stars with Mns = 0.05 in an isotropic background, at the point of closest 
binary separation along our sequence. The left and right panels show the irrotational and corotating cases, respectively. The 
position of the density maximum is -Xhs/A/bh = H-6 for the irrotational case and -Xns/Mbh = 11.8 for the corotating one. 



more detail in Section fVj As before we do not extrap- 
olate for Mns = 0.15 because the closest separation we 
can calculate is still far from the mass-shedding point. 

Since the indicator x does not exceed unity for any se- 
quence in an isotropic background and decreases mono- 
tonically toward zero, it is convenient to compare se- 
quences for different neutron star masses Mns m the 
same figure. In Fig. 1131 we show the mass-shedding in- 
dicator, Xj as a function of X-^s/r e , for several different 
sequences. Interestingly, we see that x is almost indepen- 
dent of the neutron star mass or spin, determined almost 
entirely by the normalized separation -Xns/'V A more 



massive neutron star has slightly smaller value of Xi indi- 
cating that it would be disrupted at a larger normalized 
separation ATns / r e . We can see the same behavior in the 
binary neutron stars, for example, in Fig. 7 of [38]. 



The central energy density again decreases for all se- 
quences with the binary separation, as we found for all 
sequences in a Kerr-Schild background. 
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FIG. 10: Same as Fig. [9] but for -Mns = 0.1. The position of the density maximum is X^s/Mbh = 5.24 for the irrotational 
case and Xns/Mbh = 5.37 for the corotating one. 
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FIG. 11: Same as Fig. [9] but for Mns = 0.15. The position of the density maximum is Xns/Mbh = 3.96 for the irrotational 
case and X^s/Mbh = 4.20 for the corotating one. 



V. COMPARISON BETWEEN KERR-SCHILD 
AND ISOTROPIC BACKGROUNDS 

Kerr-Schild coordinates and isotropic coordinates rep- 
resent a Schwarzschild black hole in two different coor- 
dinate systems that differ not only in the spatial coordi- 
nates but also in the time-coordinate. Spatial slices in 
the two coordinate systems therefore represent distinct 
slices of Schwarzschild, and there is no reason to expect 
that solving the conformal thin-sandwich equations on 
these different slices will result in physically equivalent 
solutions. 

It is of interest, then, to compare our results for Kerr- 



Schild and isotropic backgrounds. Most interesting, of 
course, are comparisons of coordinate-independent quan- 
tities, in particular the angular velocity (for example at 
tidal break-up), and the maximum density. We also in- 
clude comparisons of the effective enthalpy and the shift 
vector, which behaves very differently in the two coordi- 
nate systems. 
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Isotropic Isotropic Isotropic 




FIG. 12: Mass-shedding indicator x> defined by Eq. (|52p . as a function of the orbital separation for neutron star of mass 
(a)MNS = 0.05, (b)MNS = 0.1 and (c)Mns = 0.15 in isotropic backgrounds. The thick solid line with circles represents 
the irrotational case and the thick dashed line with squares the corotating one. The thin dotted and dot-dashed lines are 
extrapolations of the results for the irrotational and corotating cases, respectively. 
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FIG. 13: Mass-shedding indicator \ as a function of the orbital separation Xns/r e . Solid, dashed, and dot-dashed lines denote 
the irrotational case with neutron star masses Mns = 0.05, 0.1, and 0.15, respectively. Dotted, long-dashed, and dot-dot-dashed 
lines are for the corotating case with Mns = 0.05, 0.1, and 0.15. 



A. Orbital angular velocity at the mass-shedding 
point 

As we explained previously, we stop the computation 
of each sequence before the mass- shedding point since it 
is impossible for the spectral method to treat cusp-like 
figures. Thus, we have to extrapolate all sequences to the 
mass-shedding point to estimate the value of the orbital 
angular velocity there. In Fig. [Ml we show the mass- 
shedding indicator x as a function of the orbital angular 
velocity. The indicator \ decreases rapidly as the sep- 
aration decreases, or equivalently as the orbital angular 
velocity increases, so that extrapolating the curves using 
only our calculated values at large and medium separa- 
tions may give predictions with large errors, especially 
for the irrotational case in the Kerr-Schild background 



with Mns = 0.05; in cases where the sequences continue 
to smaller separations we expect an error of no more than 
~ 10% in our measure of the terminal angular velocity. 

For both Kerr-Schild and isotropic backgrounds, we 
find a value OMbh ~ 0.025 for the corotating case with 
neutron star mass Mns — 0.05, in good agreement with 
the results of BSS, who find OA/ B h = 0.0241 (see their 
Table 1). 

Finally, we comment on the discrepancy in the extrap- 
olation curve for the irrotational case in the Kerr-Schild 
background with Mns = 0.05. In this case, since we do 
not have results in the range < \ < 0.5, we cannot draw 
an accurate extrapolation curve, as mentioned above, and 
thus cannot use it to predict the orbital angular velocity 
at the mass-shedding point with a high degree of cer- 
tainty. However, we can predict that the real orbital 
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angular velocity at the mass-shedding point should be 
smaller than QMbh ~ 0.027 the value given by Fig. Q3J 
because the parameter \ decreases more rapidly at larger 
orbital angular velocities for the irrotational case as com- 
pared to the corotating one for the range \ < 0.5 (see 
Section V.B of 36]). This means that the data points for 
the irrotational case of the Kerr-Schild background with 
Mns = 0.05 are most likely to lie below the extrapolated 
line (solid line) in reality. 

B. Maximum value of the density quantity 

The density parameter q is a measure of the density 
as seen by a comoving observer, and is therefore gauge- 
invariant. In Fig. [IS] we show its maximum value in the 
star as a function of the orbital angular velocity, which 
is also coordinate-independent. As we have discussed be- 
fore, the maximum density decreases with binary sepa- 
ration in all cases, and decreases more rapidly for coro- 
tating configurations. 

For Mns = 0.05 we find only very small differences be- 
tween Kerr-Schild and isotropic backgrounds (less than 
1%). However, for Mns = 0.1 the compaction parame- 
ter M^ s /r (M^ s is the gravitational mass of an isolated 
neutron star and ro its radius) and hence relativistic ef- 
fects are significantly larger, the differences increase no- 
ticeably. The difference between the two backgrounds is 
now in the order of 5%, meaning that the effect of the 
background is almost as large as that of the presence of 
a binary companion. We find in all cases the Kerr-Schild 
background leads to a greater decrease in the central den- 
sity than an isotropic background. 

C. Effective enthalpy field 

For corotating configurations the definition of the en- 
thalpy can be extended to regions outside the neutron 
star by virtue of Eq. [)28p. As discussed in BSS this "ef- 
fective enthalpy field" plays the role of an effective po- 
tential, allows for the definition of a Roche lobe in a fully 
relativistic context, and is very useful for locating the on- 
set of tidal disruption. To see this, we evaluate Q28p for 
corotating binaries 

flcnt = C-i/ + ln 7o , (54) 

where C is a constant. Inserting the definition of v (Eq. 
©) and 70 (Eq. (f2"4"| ) for the corotating case, we obtain 

ff ent = C-lna-iln(l-^^). (55) 

In the Newtonian limit this reduces to 

H^ = C'-d> + ^(Qxr) 2 , (56) 

where <j> denotes the total gravitational field and C is an 
integration constant which is the Newtonian limit of C. 



For irrotational configurations Eq. (|28[) depends on the 
solution of the continuity equation, which can only be 
solved in the stellar interior. In this case we cannot ex- 
tend the definition of the enthalpy to regions outside the 
neutron star, meaning that a straight-forward definition 
of the effective potential is possible only for corotating 
configurations. 

We show contours of the effective enthalpy field (j55"|) 
in Figs. [TBI and [T7l for neutron star of mass Mns = 0.05 
and 0.1, respectively. In each figure, the left panel is for 
the Kerr-Schild background and the right panel for the 
isotropic one. In both figures, we show the position of 
the inner Lagrange point using the symbol "x" and that 
of the outer one with "+". The equipotential surface 
passing through the inner Lagrange point defines the rel- 
ativistic Roche lobe. It is clear that the neutron star in 
our results still fits well within its Roche lobe, but if we 
decrease the separation only a few percent further, we 
expect a rapid deformation on the inner edge as a cusp 
forms at the mass-shedding limit. 

For the Kerr-Schild background the absence of a sym- 
metry across the X—Z plane is again quite noticeable, in 
that the Lagrange points do not lie on the X-axis. For 
an isotropic background, on the other hand, the presence 
of this symmetry force the Lagrange points to lie on the 
A- axis. 



D. Shift vector 

One of the most significant differences between the 
Kerr-Schild coordinate system and the isotropic one is 
the existence/ absence of the black hole shift vector. 
In isotropic coordinates the background shift is zero, 
while for a Kerr-Schild background the main contribu- 
tion to the total shift arises from the Schwarzschild back- 
ground. We show the shift vector for an inertial observer, 
/^ns+^bh, for corotating configurations with Mns = 0.05 
in Fig. [THJ The left panel is for the Kerr-Schild case at 
the orbital separation Ans/Mbh = 13.0 and the right 
panel for the isotropic case at Ans/Mbh = H-8. The 
neutron star rotates counterclockwise. It is evident that 
for the Kerr-Schild background the shift is dominated by 
the outward-pointing background contribution /3 BH . The 
lack of a symmetry across the X—Z plane is again quite 
obvious for the Kerr-Schild background. 



VI. SUMMARY 

We have computed quasiequilibrium sequences of 
BHNS binary systems in general relativity. Under the 
assumption of an extreme mass ratio, Mbh 3* Mms, 
we have treated the contribution of the black hole 
gravitational field as a fixed background, adopting the 
Schwarzschild metric in both Kerr-Schild and isotropic 
coordinates. The neutron star is modeled as a relativis- 
tic polytropic EOS with adiabatic index T = 2. We have 
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FIG. 14: The mass-shedding indicator \ as a function of the orbital angular velocity. The left panel shows sequences with 
neutron star mass Mns = 0.05, the right Mns = 0.1. K-S and ISO denote sequences performed in the respective backgrounds. 




FIG. 15: Maximum value of the density quantity g max as a function of the orbital angular velocity. The left panel corresponds 
to a neutron star of mass Mns = 0.05 the right panel A/ns = 0.1. K-S and ISO denote the respective backgrounds. 



solved a set of equations for two rotation states, irro- 
tation and corotation. These results generalize those of 
BSS, which presented the first relativistic, self-consistent 
method to evaluate quasiequilibrium black hole-neutron 
star configurations. 

We have developed a new numerical code for the 
present study, based on a spectral methods numerical 
approach that was used for a series of works on binary 
neutron stars [H], [H, H3, [H, H^, H(| • After confirming the 
validity of the numerical code in several tests, we have 
computed quasiequilibrium sequences of the BHNS bina- 
ries with a mass ratio Mbk/Mns = 10 and neutron star 
masses Mns = 0.05, 0.1, and 0.15. We expect that spec- 
tral methods allow us to construct each individual model 
with greater accuracy than the finite difference methods 



used by BSS. However, the assumption of smoothness 
inherent in spectral methods prevents us from construct- 
ing models very close to the onset of tidal disruption. 
This leads to errors in parameters describing the onset 
of tidal disruption that are greater than those obtained 
with finite differencing. While we maintain here some of 
the assumptions from BSS, including extreme mass ra- 
tios and polytropic equations of state, and plan to relax 
these in the near future. 

In agreement with earlier studies (e.g. [27[), we find 
that the effect of rotation on the onset of tidal disruption 
is fairly small and in the order of a few percent. This does 
not rule out, however, that rotation has a larger effect 
on the dynamics of the tidal disruption itself (compare 
[IH). In particular, the rotational state of the neutron 
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FIG. 16: Contours of the enthalpy field H Bnt extended outside of the neutron star, which has a mass Mns = 0.05 and is 
corotating. The thick solid curve denotes the stellar surface, thin solid curves are contours located inside the neutron star, 
and dashed curves those outside. The left and right panels show configurations in the Kerr-Schild and isotropic backgrounds, 
respectively. Both depict the innermost point along the respective sequences, listed in the last line of the corresponding tables. 
The symbols x and + denote the inner and outer Lagrange points, respectively. Note that the horizontal axis X in the left 
panel is the X coordinate of the Kerr-Schild coordinate system and that in the right panel is in the isotropic coordinate system. 



star may well affect the size (or existence) of an accretion 
disk that may form as the neutron star is disrupted. Such 
an accretion disk is at the core of BHNS models as central 
engines of short-period gamma ray bursts (e.g. [1, H, [l(| ) . 

We also find that the choice of background does have 
some effect on coordinate- independent quantities describ- 
ing the resulting binary configurations. This indicates 
that the choice of background may affect the degree 
to which the solutions approximate quasiequilibrium, as 
well as the amount of spurious gravitational radiation 
inherent in the solutions. It will be very interesting to 
study these differences in future dynamical simulations 
(compare (l4|)- 
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APPENDIX A: EQUATIONS 



In this appendix, we present the explicit forms of equations which we solve in our numerical code. Having applied 
the decompositions of metric quantities (fT5)) . ([TO]) , and (j2"D|) to the equations {7}, ([3]), and (fTU|) , we can derive the 
final forms for the neutron star components in a numerically convenient form, inserting the black hole components as 
a fixed background. 



Kerr— Schild (Corotation) 



Isotropic (Corotation) 




FIG. 17: Same as Fig. 1161 but for a neutron star of mass Mns = 0.1. 



1. Kerr-Schild background 



a. Gravitational field equations 



The equations (J7J), (g]), and (fTU|) , respectively, can be written as 
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Kerr— Schild (Corotation) 



Isotropic (Corotation) 
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FIG. 18: Shift vectors seen by an inertial observer in the X—Y plane (equatorial plane). The left and right panels show the 
Kerr-Schild and isotropic backgrounds, respectively. In both cases the neutron star has a mass M^s = 0.05, and is corotating. 
Thick solid circles in each panel denote the stellar surface, and the arrows the direction of the shift vector. Note that the 
horizontal axis X in the left panel is the X coordinate of the Kerr-Schild coordinate system and that in the right panel is in 
the isotropic coordinate system. 
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where A and V» are the fiat Laplace operator and partial derivative. The conformally related trace-free extrinsic 
curvature of the neutron star part is defined as 
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b. Equation of continuity 

Having denned a velocity field related to the orbital motion, i.e., 

W* = ^h lu Ul 

we can express the velocity potential as 

* = * +%W V, 
where Wq is the constant value of W l at the center of the neutron star, i.e., 

Wq = (tp 4 hj n Uo) con tor = constant. 



(A5) 



(A6) 



(A7) 



Since the Newtonian limit of the right hand side of Eq. (|A5|) is the orbital motion (f2 x tbh) 1 , Eq. (|A7[) is merely the 
relativistic analogue of the translational motion of the center of the neutron star. Thus, the quantity ^g is regarded 
as the residual of the velocity potential, once the constant rotational component is subtracted away. The gradient of 
\l/g yields the counter-rotation seen by a co-orbiting observer. The gradient and the Laplacian of vPq become 



V 1 * = V 4 * + ^g, 

A* = Afg. 



(A8) 
(A9) 



Inserting the Kerr-Schild metric for the black hole into Eq. 



we have our final form of the equation of continuity, 
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c. Determination of the orbital angular velocity 
The left-hand side of the equation P2"|) can be expressed 
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(All) 



where (3[ denotes the shift vector seen by the inertial observer, i.e., j3\ = /J^s + /3 BH . 

The orbital angular velocity appears explicitly in the above equation, which corresponds to the left-hand side of 
Eq. (|42[) . On the contrary, we compute the right-hand side of Eq. (|42p in the form shown in which it is shown. 
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2. Isotropic background 



a. Gravitational field equations 



In the case of isotropic backgrounds, Rij = R = Ag H = K = 0, because /3 BH = and 7^ = Tjij. This reduces the 
field equations to a much more simple form than those in Kerr-Schild backgrounds. 
The equations 0, (HJ), and (fT0|) are re- written as follows. 
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b. Equation of continuity 



The equation of continuity (I32[) is written 
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Here we have used the same quantities <£>0j W l , and Wq 1 as defined in the Kerr-Schild part, Sec. IA 1 bl 



c. Determination of the orbital angular velocity 
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Note here that since there is no black hole shift vector in the isotropic background, we write the shift vector seen by 
the inertial observer as (3^ s . 
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